On the width and shape of gaps in protoplanetary disks 
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Abstract 



Although it is well known that a massive planet opens 
a gap in a proto-planetary gaseous disk, there is no 
analytic description of the surface density profile in 
and near the gap. The simplest approach, which is 
based upon the balance between the torques due to 
the viscosity and the gravity of the planet and as- 
sumes local damping, leads to gaps with overesti- 
mated width, especially at low viscosity. Here, we 
take into account the fraction of the gravity torque 
that is evacuated by pressure supported waves. With 
a novel approach, which consists of following the fluid 
elements along their trajectories, we show that the 
flux of angular momentum carried by the waves corre- 
sponds to a pressure torque. The equilibrium profile 
of the disk is then set by the balance between grav- 
ity, viscous and pressure torques. We check that this 
balance is satisfied in numerical simulations, with a 
planet on a fixed circular orbit. We then use a ref- 
erence numerical simulation to get an ansatz for the 
pressure torque, that yields gap profiles for any value 
of the disk viscosity, pressure scale height and planet 
to primary mass ratio. Those are in good agreement 
with profiles obtained in numerical simulations over 



a wide range of parameters. Finally, we provide a 
gap opening criterion that simultaneously involves 
the planet mass, the disk viscosity and the aspect 
ratio. 

1 Introduction 

The dynamical evolution of planets in proto- 
planetary disks has become a topic of renewed 
interest in the last decade, boosted by the discov- 
ery of extra-solar planets, and in particular of hot 
Jupiters. In fact, the observation of giant planets 
close to their parent stars argues for the existence 
of effective mechanisms of planetary migration, 
which can be found in the study of planet-disk 
interactions. 

Several types of migration have been identi- 
fied, depending on how the planet modifies the 
local density of the disk. Type I migration oc- 
curs when the planet is not massive enough to 
significantly alter the local density of the disk; 
the planet migrates inward with a speed propor- 
tional to its mass (Ward, 1997). Type II migra- 
tion corresponds to the case where the planet is so 
massive that it opens a clear gap in the disk ; the 
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migration then depends on the viscous evolution 
of the disk (Lin and Papaloizou, 1986 a, b). Type 
III (or runaway) migration corresponds to plan- 
ets with intermediate mass, which do not open 
a clear gap, but only form a dip around their or- 
bits in the gas surface density profile ; under some 
conditions, their migration drift rate can grow ex- 
ponentially, in a runaway process (Masset and Pa- 
paloizou, 2003). 

The modification of the disk density is the re- 
sult of the competition of torques exerted on the 
disk by the planet and by the disk itself. More 
precisely, the planet gives some angular momen- 
tum to the outer part of the disk, while it takes 
some from the inner part (Lin and Papaloizou, 
1979; Goldreich and Tremaine, 1980). In doing 
so, it pushes the outer part of the disk outward 
and the inner part inward, and therefore tends to 
open a gap. However, the internal evolution of 
the disk, which tends to spread the gas into the 
void regions, opposes to the opening of the gap. 

However, there is a lack of an analytical predic- 
tion of the gap profiles. Classically, the gap is con- 
sidered to have a step function profile, with the 
edges located at the sites where the total gravity 
torque is equal to the total viscous torque. This 
is obviously an oversimplification. A more sophis- 
ticated approach has been recently presented by 
Varniere et al. (2004). They provide an ana- 
lytic expression that describes the gap profile, by 
equating the viscous and gravity torques on any 
elementary disk ring. We will provide more de- 
tails on this approach in section |21 The problem 
is that, when the viscosity is small, the viscous 
torque is small as well, and thus it cannot coun- 
terbalance the gravity torque. Consequently, in a 
low viscosity disk, a non-migrating planet should 
open a very wide gap, unlike what is observed in 
numerical simulations (see e.g. Fig. [TJ : gaps do 
increase in width and depth as the viscosity de- 
creases, but the dependence of the gap profile is 
less sensitive on viscosity in the numerical simu- 
lations than it is expected in theory. 

The reason for this difference is that not all 
of the gravity torque is locally deposited in the 
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Figure 1: Gap profiles created by a Jupiter mass 
planet, for different viscosities. The vertical axis 
represents the normalized azimuth-averaged density. 
The horizontal axis represents the distance to the 
primary in normalized units. Top panel : analytic 
curves obtained by matching the differential torques 
due to gravity and viscosity on each elementary 
ring of the disk. Bottom panel : numerical profiles 
obtained in simulations, after 1000 planetary orbits 
for the three largest viscosities, and 5000 orbits for 
log(i^) = —5.5 and —6.5. 

disk. It is transported away by density waves 
(Goldreich and Nicholson, 1989 ; Papaloizou and 
Lin, 1984 ; Rafikov, 2002 ; see Appendix C). These 
waves are observed in simulations. In this situ- 
ation, the viscous torque has to counterbalance 
only a fraction of the gravity torque, which yields 
narrower gaps than expected from the simple vis- 
cous/gravitational torque balance. 

An evaluation of the fraction of the gravita- 
tional torque that is locahy deposited at shock 
sites in the disk has been undertaken by Rafikov 
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(2002). However he did not use his analysis to 
provide an analytic representation of the gap pro- 
file. Moreover, his calculation required several as- 
sumptions (the planet Hill radius had to be much 
smaller than the local disk height, the disk surface 
density was assumed to be uniform, etc.), which 
are not satisfied in the general giant planet case. 

Here we introduce a novel approach. We fol- 
low a fluid element along its trajectory which, in 
the steady state, is periodic in the planet corotat- 
ing frame. A flux of angular momentum carried 
by the density waves corresponds to a pressure 
torque acting on the fluid element, whose average 
over a synodic period is non-zero. In this work, we 
evaluate this averaged pressure torque, together 
with the gravity and viscous torques. The fact 
that the fluid element path is closed implies that 
these time averaged torques balance. 

In section 0] we introduce the pressure torque, 
and we use numerical simulations to check that 
the gap structure is set by the equilibrium be- 
tween the gravity torque and the sum of the vis- 
cous and the pressure torques. In section El we 
construct a semi-analytic algorithm and we get 
an expression to compute that gap profile. We 
compare our results with the profiles obtained in 
numerical simulation, and we discuss the merits 
and limitations of our method. In section El we 
use our algorithm to explore the dependence of 
the gap structure on disk viscosity and aspect ra- 
tio. We recover the trends observed in numer- 
ical simulations, namely the limited gap width 
in low viscosity disks and the filling of the gap 
with increasing viscosity and/or aspect ratio. Fi- 
nally in section ini we provide a gap opening cri- 
terion that simultaneously involves the viscosity, 
the scale height and the planet mass. 

2 Gravity and Viscous 
Torques 

In this section we revisit the calculation of the 
gravity and viscous torques mentioned in the 
Introduction. We show that considering them 



alone, as usually done, is not sufficient to achieve 
a quantitatively correct description of the gap 
profiles. 

2.1 Notations 

The disk is represented in cylindrical coordinates 
{r,6,z), centered on the star, where the plane 
{z = 0} corresponds to the mid-plane of the disk. 
The disk viscosity z/ and aspect ratio (H/r) - 
where H denotes the thickness of the disk- are 
assumed to be invariant in time and space. The 
equations of fluid dynamics are integrated with 
respect to the z-coordinate, so that z disappears 
from the equations and only two dimensions are 
effectively used. This procedure introduces the 
concept of surface density S, which is defined as 
J^l^ pdz, where p is the volume density in the 
disk. 

In the theoretical analysis (but not in the nu- 
merical calculations) the disk is assumed to be 
axisymmetric, so that E only depends on r. The 
angular velocity Q is assumed to be Keplerian : 
n oc r~^/^. The planet is assumed on a circular 
orbit around the star. The radius of its orbit is 
denoted r^. The mass of the planet is denoted 
Mp and its ratio with the mass of the central star 
M* is q. Normalized units are introduced, so that 
M* = Vp = 1 and the gravitational constant G 
is also assumed to be unity. In the limit g — >■ 0, 
this sets the angular orbital velocity of the planet 
Qp = 1 and its period equal to 27r. 

2.2 Total torques 

Usually, one considers the part of the disk extend- 
ing from a given radius ro > Vp to infinity. The 
study of the part of the disk extending from to 
ro < rp is done in an analogous way. Two torques 
are evaluated. The first one is due to the disk vis- 
cosity and can be easily derived from the stress 
tensor in a Keplerian disk with circular orbits. 
The torque exerted on the considered part of the 
disk (r > ro) by the complementary part is writ- 
ten (see for instance Lin and Papaloizou, 1993) : 
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T,^ = SnTiuro'^Qo (1) 

Notice that more refined expression for perturbed 
disks with eccentric orbits have been proposed in 
the hterature (see for instance Borderies et al. 
1982), but they have not been used in the works 
that we review in this section. 

The second torque comes from the gravity of 
the planet. It can be computed following two 
different approaches. In the first one (Goldreich 
and Tremaine, 1980; Ward, 1986), it is decom- 
posed into the sum of the individual torques ex- 
erted at each Lindblad resonance. In the second 
approach (Lin and Papaloizou, 1979 ; Goldreich 
and Tremaine, 1980) it is obtained by computing 
the angular momentum change for fluid elements 
at conjunction with the planet, using an impulse 
approximation. The two approaches are known to 
give equivalent results. In the following, we use 
the expression from the impulse approximation : 



Thus, a gap can be opened only if : 



Ar 



(2) 



where Aq = {j-q — r-p). The above expression 
gives the torque exerted by the planet on a disk 
extended from tq to inflnity. It is valid only 
for |Ao| > Am, where Am is the maximum of 
H (the local thickness of the disk) and the Hill 
radius of the planet Rh = {q/^Y^^ (Goldreich 
and Tremaine, 1980; Ward, 1997). The value of 
the numerical coefficient C depends on the ap- 
proach followed for the calculation of the torque. 
In the most recent and refined calculation, Lin 
and Papaloizou (1993) found C = ^[2ifo(i) + 
-^i(i)]^ ~ 0.836 (where Kq and Ki are modified 
Bessel functions). 

Classically, the gap is modeled step func- 
tion profile in the disk surface density, with edges 
placed at a distance Aq from the planet orbit, 
with Aq given by the solution of the equation 
Tg(Ao) = T,^, and Tg and given in ^ and 
(PJ). The maximal gravity torque is: 



Tg(Am) ~ 0.836 g^Sr/fij 



A. 



u < 0.0887g 



2 ' p ^''p 



(3) 



otherwise T^, is larger than Tg and the gas over- 
runs the planet. Condition (jH} is equivalent to 
that given in Bryden et al. (1999), expressed as 
a constraint on the mass of the planet relative to 
the viscosity of the disk. 

2.3 Differential torques and com- 
parison with numerical tests 

Differential torques. Varniere et al. (2004) 
proposed a more refined approach to model the 
surface density profile in the gap. Their approach 
is based on a simple consideration : in equilib- 
rium, when a steady state is reached, the gravity 
torque and the viscous torque must be equal on 
every elementary ring of the disk. The torques 
acting on elementary rings can be computed by 
differentiation relative to r = tq of (0) and Q : 
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r dS 1' 
S d7 ^ 2 



(27rrS) (4) 



STg{r) ^ 0.4 gV/f]/r-i (^^) ^ (27rrS) (5) 



Matching 5Ty and 6Tg gives a differential equation 
in S : 

1 dE ^ 6Tg{r) _ 1_ 

The integration of this equation gives the profile 
of the gap. 

The top panel of Fig. ^ gives examples of the 
solution of Eq. (jH)) for several values of the vis- 
cosity, from strong (z/ = 10~^'^^) to weak (z/ = 
10~^'^). To compute them from (jU)), we have (i) 
assumed that the mass of the planet is 10~^ in our 
normalized units, (ii) imposed the boundary con- 
dition S(ro = 3) = l/y/S and (iii) assumed that 
the gravity torque is null in the horseshoe region, 
here approximated by : rp—2RH <r < rp+2RH. 
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As a consequence of (iii) the surface density pro- 
file in the vicinity of the planet assumes an equi- 
librium slope proportional to 1/ y^, which makes 
STy null. Notice that the slopes of the surface 
density at the edges of the gap do not depend 
on our assumptions (ii) and (iii), but are dictated 
solely by the differential equation We re- 

mark that the profiles illustrated in the figure are 
the same as in Varniere et al. 2004, despite the 
fact that these authors consider the gravity torque 
as given by the sum of the individual Lindblad 
resonances. This again underlines the equivalence 
of the two approaches for the calculation of the 
gravity torque discussed in 12.21 

Numerical simulations. We have tested the 
results of these analytic calculations using purely 
numerical simulations. For this purpose, we have 
used the 2D hydrodynamic code described in 
Masset (2000), and considered a Jupiter mass 
planet (g = 10~^) in a disk, whose initial sur- 
face density profile decays as I/a/t, and S(rp) = 
6.10"^ (the value of the minimal mass solar neb- 
ula at 5 AU (Hayashi, 1981)). The disk aspect 
ratio was fixed at 5%. The viscosity was chosen 
equal to the values used for the analytic computa- 
tions, for direct comparison. In these simulations, 
the planet was assumed not to feel the gravity of 
the disk, so that it did not migrate. The grid used 
by the code for the hydrodynamical calculations 
extended from 0.5 to 3 (we remind that the planet 
location is Vp = 1). The boundary conditions in 
r are non-reflecting, which means that the waves 
behave as if they were propagating outside the 
boundaries of the grid. The angular momentum 
they carry is thus lost ; we have checked that the 
flux through the outer boundary represents only a 
neghgible fraction of the total gravity torque (see 
Appendix C). The relative surface density ampli- 
tude perturbation at the outer boundary has been 
measured to be less than 5%. The size of the grid 
was 150 cells in radius and 325 cells in azimuth. 
The simulations were carried on for 1000 plane- 
tary orbits, for viscosity down to 10~^'° and 5000 
orbits for weaker viscosities. At these times, the 



profile of the gap does not seem to evolve signif- 
icantly any more, as also found by Varniere et al. 
(2004). 



Comparisons. The results are illustrated in 
the bottom panel of Fig. ^ As anticipated in 
the introduction, we remark an evident differ- 
ence with the analytic predictions. The simu- 
lated gap is much narrower than the one predicted 
by the analytic expression © for low viscosities 
{v < 10-6). At first sight, one might think that 
the discrepancy between the analytical and nu- 
merical solutions is due to the numerical viscos- 
ity (dissipation due to numerical errors) of the 
computer code. However, this is unlikely for the 
following reasons : (i) different gap profiles are ob- 
served for different viscosities, which shows that 
the simulation is not dominated by the numer- 
ical viscosity, as the latter should be the same 
in all simulations ; (ii) changing the resolution 
of the grid used in the numerical scheme, which 
changes the numerical viscosity, does not affect 
the gap profiles significantly; (iii) different nu- 
merical schemes give consistent results (De Val- 
borro, private communication). 

As anticipated in the introduction, the problem 
with this analytical modeling is the assumption 
that the gravity torque is entirely deposited in 
each annulus of the disk. A condition for such de- 
position to happen is that > H (Lin and Pa- 
paloizou, 1993). Thus, this is usually considered 
second independent criterion for gap opening, 
in addition to Q (Bate et al, 2003). However, 
even if this condition is satisfied, a fraction of 
the gravity torque is still evacuated by the waves 
(Goldreich and Nicholson, 1989 ; Papaloizou and 
Lin, 1984 ; Rafikov, 2002). The problem is to eval- 
uate this fraction. Below we show that it can be 
computed from a mean pressure torque acting on 
the fluid elements over their periodic equilibrium 
trajectories. 
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3 Pressure torque 

Consider an arbitrary closed curve in the disk and 
a little tube around it. The rate of change of an- 
gular momentum of the matter in the tube is the 
sum of the differential flux of angular momentum 
through its boundaries (due to the advection of 
matter) and of the torques acting on it. From the 
Navier-Stokes equations, in addition to the grav- 
ity and viscous torques, there is a third torque 
due to pressure : 

tp = fc,'^d9, (7) 

where the integral is computed along the curve 
and we have assumed the usual state equation 
P = Cs^S (with Cs = HQ denoting the sound 
speed). If the tube is a ring centered at the origin 
r = 0, the torque tp is equal to zero (on the ring, 
r is constant and dTi/dO = d'E/dO), while the dif- 
ferential flux of angular momentum is generally 
not zero. The latter is the flux carried by the 
pressure supported wave (Goldreich and Nichol- 
son, 1989; see Appendix C). On the contrary, if 
one chooses a stream tube {i.e. a tube bounded 
by two neighboring streamlines), the differential 
flux is obviously zero (there is no flux of mat- 
ter, by definition of stream tube), while tp is non 
zero in general. The latter is true because the 
streamlines are strongly distorted (see Fig. ^ so 
that S = E{r{9),9) and thus dE/dO ^ dS/d^. 
This shows that one can translate the angular mo- 
mentum flux carried by the waves into a pressure 
torque, by a suitable partition of the disk in con- 
centric tubes. Obviously the two approaches are 
equivalent as the physics is the same. However, 
working with stream tubes and pressure torques 
gives practical computational advantages. This 
is therefore the approach that we follow in this 
paper. 

Below, we check that the gravity, pressure and 
viscous torques really cancel each other in the 
numerical simulations, once the steady state is 
reached. 




Figure 2: Disk surface density map in the vicinity 
of a gap opened by a Jupiter mass planet located at 
(rp = l,9p = 0). Light grey denotes high density and 
black low density, in a logarithmic scale. The white 
curves show some streamlines, in the frame corotat- 
ing with the planet. They are followed from vr to — vr 
for r > 1, and from — vr to vr for r < 1, periodically. 
Two of them correspond to horseshoe orbits in the 
planet corotation region. Notice the strong distortion 
of the streamlines when they cross the over-density 
corresponding to the spiral wave (wake) launched by 
the planet. 

3.1 Computation of the torques 
along the trajectories 

The approach outlined above requires that the 
stream tubes are closed. In our simulations this 
is true at the steady state, because our boundary 
conditions preserve the initial radial velocity at 
the edges of the grid, which is null result of 
our choice for the initial disk density profile S oc 
1 / ^/r. The calculation of the streamlines, which 
is done in Fourier space to ensure periodicity, is 
detailed in Appendix A. We remark that, in the 
steady state, the streamlines coincide with the 
fluid element trajectories. 

Denoting the streamline by rj(6'), we numeri- 
cally compute the following expressions, which 
are the integrals of (l/S)(rFe) with Fg the az- 
imuthal component of the force due to gravity. 
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viscosity, or pressure respectively : 



n{e'),e') 



09' 



dO' 



9' 
(9) 



tp{9) 



W) 



9S 



89' 



d9' 



(10) 

Here, denotes the gravitational potential of the 
planet, and = 



where T 



T . 



is the local viscous 



stress tensor for a Newtonian fluid : T = 
2Sz/ {j) — (|V'y)/j, where D is the strain tensor 

and / is the identity matrix. We integrate from 
TT to 9, with > 9 ^ — vr, because we consider 
To > Tp, so that the angular velocity is negative in 
the corotating frame. The time required to reach 
9 from TT is denoted Ti{9). Thus, as the trajec- 
tories coincide with the streamlines in the steady 
state, the expressions above describe the averaged 
torques felt by a fluid element that travels from 
the planet opposition to 9. 

In the following, we denote for simplicity by 
tg, tu, tp the expressions (jH|)- (fTIH) evaluated at 
9 = — TT. The total torques acting on the stream 
tube centered around the considered streamline 
are simply the product of these quantities times 
the mass carried by the tube. 

In the next paragraphs we give a brief descrip- 
tion of the integrated torques (jH|)- (fTm) as functions 
of 9, which are plotted in Fig.Elfor the streamline 
starting at r = 1.2 at opposition with the planet. 

viscous torque : The growth of the integrated 
viscous torque appears to be nearly linear with 
respect to the azimuth, leading to a total nega- 
tive torque. We verified that on this streamline 
^ 5T^/(27rSr), with 6T^ from 0. Thus, the 
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Figure 3: Graphical representation of the expres- 
sions © (integrated gravity torque, bold short-dash 
curve), (jU (integrated viscous torque, bold dotted 
curve) and (|1U|) (integrated pressure torque, bold 
solid curve). Their reference scale is reported on the 
left vertical axis. The streamline followed for their 
calculation is plotted in the planet corotating frame 
dashed curve at the top of the figure and can 
also be seen on Fig. |2 while the position of the planet 
is shown by a filled dot at the bottom ; the corre- 
sponding scale is reported on the right vertical axis. 
The thin solid curve shows the difference between the 
angular momentum measured along the streamline 
{H{6) ), and the sum of the three integrated torques 
and of the initial angular momentum {7i{6) ). A 
small difference is almost impulsively acquired at the 
wake crossing, due to numerical approximations. 



viscous torque depends only on the radial rela- 
tive derivative of the azimuthally averaged den- 
sity However, on streamlines that pass 
closer to the planet, the difference between ti, and 
5Tjj/ (27rSr) becomes more significant (see Fig.lHl). 

gravity torque : The evolution of this inte- 
grated torque is not monotonic. The fluid ele- 
ment is first repelled by the planet, as a result of 
the indirect term in the gravitational potential. 
Then, when 9 decreases below ~ 0.5, it starts to 
be attracted by the planet. The attraction be- 
comes stronger and stronger as the fluid element 
approaches conjunction, namely as 9 decreases to 
0. The integrated torque becomes negative. After 
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conjunction, the planet tends to pull the fluid el- 
ement toward positive 9, giving a positive torque. 
As a result, the fluid element is rapidly repelled 
toward larger r, as one can see from the tra- 
jectory on Fig. ISl This is typical of the scatter- 
ing of test particles in the restricted three body 
problem, which quahtatively justifles the impulse 
approach for the calculation of the gravitational 
torque, as in Lin and Papaloizou (1979). 

However, Lin and Papaloizou's calculation 
holds in the approximation r ~ rp. By compar- 
ing the numerical estimate of 5Tg/(27rSr) with 
5Tg given by Eq. (0), we flnd that the follow- 
ing expression, which has the same dependence 
in A and nearly the same numerical coefficient, 
but which distinguishes r and r^, provides a 
much more accurate representation of the grav- 
ity torque : 

= 0.35 gV/ f]/ r sgn(A) , (11) 

In reality, tg depends on the exact shape of the 
streamlines, which in turn depends on the scale 
height and the viscosity (see Fig. HI and EI). How- 
ever, the difference is moderate and limited to the 
vicinity of the planet, so that in the following we 
use expression (fTT|) for all cases. 

We stress that expression (fTTj) gives the torque 
exerted on the fluid element, which is generally 
not the torque deposited in the disk. In fact, even 
in the absence of viscosity, it does not correspond 
to the change of angular momentum of the fluid 
element, because some of the angular momentum 
is carried away by the pressure torque. 

pressure torque : The dependence of this 
torque on 9 is simple to understand if one takes 
into account that : (i) the trajectories cross the 
wake immediately after the conjunction with the 
planet (Fig. |21) ; (ii) the wake is a strong over- 
density in the disk ; (iii) the pressure term 9S / 89 
makes over-densities repellent. Thus, as the fluid 
element approaches the wake, its azimuth 9 de- 
creasing in the corotating frame, the pressure rises 
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Figure 4: Absolute value of tg as a function of r, 
where r is the distance to the central star at planet 
opposition of the streamline on which the gravity 
torque is integrated. The bold long-dashed curve 
is obtained from the simulation with q = 10~^, 
{H/r) = 5%, and u = 10~^-^, while the bold short- 
dashed curve corresponds to a more viscous case 
[v = 10~^'^). The solid line traces expression (|11() . 
which remarkably fits the results obtained in the less 
viscous case. The dashed-dotted curve traces expres- 
sion which shows that 5Tg ^ tg. 

and tends to push the fluid element back in the di- 
rection of increasing 9. This gives a positive local 
torque and it explains the peak in the integrated 
pressure torque in Fig. El Then, after that the 
fluid element has crossed the wake, the pressure 
decreases as 9 decreases. This leads to a negative 
local pressure torque. It corresponds to the fall 
after the peak on Fig. [HI The negative contribu- 
tion is bigger than the positive one because of the 
asymmetry of the trajectory relative to the wake 
position, which is clearly visible in Fig. El 

Clearly, the pressure torque must depend on 
the shape of the streamlines and on the surface 
density relative radial gradient, which govern the 
shape of the wake and its density enhancement. 
We return to this in section |31 

3.2 Torque balance at equilibrium 

From Fig. |21we remark that, at = — vr, the sum 
of the viscous and pressure torques is basically 
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Figure 5: Sketch on the origin of the pressure torque. 
Here are drawn some streamlines, the width of which 
represents the mass carried by the corresponding 
streamtube. At the gap edge, the pressure gradient 
gives a force. The distortion of the streamlines at 
the wake leads to a large azimuthal component of 
this force, which gives a torque. 

the opposite of the gravity torque. Therefore, the 
three torques approximately balance out. 

Given the angular momentum H of a fluid ele- 
ment at 6 = IT, one can compute the angular mo- 
mentum Ti.{9) that it would have if its trajectory 
were governed exclusively by the three torques 
mentioned above : n{6) = H{7r) + tg{e) + U{0) + 
tp{6). This can be compared with the local angu- 
lar momentum on the trajectory H{6), measured 
directly from the numerical simulation. In Fig. El 
the thin line shows H{6) — T-i{9). This function is 
zero for 9 evolving from tt down to ~ 0, where the 
wake is crossed. At the wake crossing, a small kick 
is observed. Then, when 9 evolves from the wake 
location to — vr, the function H{9) — T-C{9) remains 
constant again. This confirms that the trajectory 
is essentially governed by the three torques men- 
tioned above. 

The small difference between H and 7i (Fig. ^ 
could in principle be due to the pseudo-viscous 



pressure introduced in the simulation to avoid nu- 
merical instabilities (Lin and Papaloizou, 1986a), 
but we have verified that the effect of the latter 
is negligible. Thus, we conclude that it is a con- 
sequence of numerical errors, introduced by the 
grid discretization at the shock site. This numer- 
ical issue evidently prevents the three cumulative 
torques from balancing out perfectly at 6* = — vr : 
indeed, their sum is equal to ?i(— vr)— i7(— tt) 7^ 0. 

In order to explore the relative importance of 
viscosity and pressure in different situations, we 
show in Fig. IHl the three averaged torques as a 
function of r for two simulations, with u = 10~^'^ 
(top panel) and u = 10~^'^ (bottom panel). In the 
more viscous case, the pressure torque becomes 
relevant for r < 1.2, i.e. at the edge of the gap. 
There, it substantially helps the viscous torque 
in counterbalancing the gravity torque. This ex- 
plains why the gap observed in the simulation 
is narrower than the one predicted by the the- 
ory considering only the gravity and the viscous 
torques alone (see Fig. IH). In fact, if the pres- 
sure torque were not present, all over the region 
r < 1.2 the relative radial gradient of the surface 
density of the disk would have needed to be much 
steeper, in order to enhance the viscous torque up 
to the value of the gravity torque (see Eq. (jU). 
This would have given a wider and deeper gap 
profile. 

It is interesting to compare the top panel of 
Fig. ini with the lower panel, which is plotted for 
a value of the viscosity that is almost an order 
of magnitude smaller. First, we remark that the 
gravity torque is somewhat smaller in the vicinity 
of the planet ; this is due to a (moderate) change 
of the shape of the streamlines, as discussed in 
last subsection. The viscous torque has decreased 
much more than the gravity torque, but not pro- 
portionally to the viscosity; this is because the 
profile of the gap has changed and the relative ra- 
dial gradient of the surface density is now steeper. 
The pressure torque has increased relative to the 
gravity torque, and is now non-negligible in the 
full region r < 1.3. It is always larger in absolute 
value than the viscous torque. Its radial profile 
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Figure 6: The torques tg, ty, tp are plotted in bold 
lines function of r, which denotes here and in 
the following plots the radius of the streamline at 
opposition with respect to the planet {9 = ibvr) ; the 
horseshoe region r < 1.15 ~ rp + 2 Rh (see Fig. ^ is 
excluded. The thin dotted line shows the value of the 
viscous torque given by 6Ti^ / {2TTTir) , with 6Ti, from 
(jlj) (keplerian circular approximation) . The thin solid 
line is the sum of the three torques. It is not exactly 
zero, in particular in the vicinity of the planet, be- 
cause of numerical approximations generated at the 
wake crossing. Top panel : large viscosity case ; the 
pressure torque becomes relevant only close to the 
planet. Bottom panel : low viscosity case ; the pres- 
sure torque appears further from planet, compensat- 
ing for the smaller viscous torque. 

looks very similar to that of the gravity torque. 
In essence, it is the pressure torque that coun- 
terbalances the action of the planet, with the vis- 
cosity only playing a minor role. Thus there is a 
dramatic qualitative change, with respect to the 
previous case, in how the torques balance out to 



settle the equilibrium configuration. 

The two cases discussed above convincingly 
show that the disk equilibrium is set by the equa- 
tion 

tg+U+tp = 0. (12) 

When the viscosity fades, the role of pressure 
takes over in controlling the gap opening process, 
limiting the gap width. This means that, as vis- 
cosity decreases, a larger fraction of the gravity 
torque is transported away by the pressure sup- 
ported waves. This phenomenon explains why 
the width of the gap increases with decreasing 
viscosity in a much less pronounced way than in 
Varniere et a/.'s model, which does not include a 
pressure torque. 

The role of pressure in limiting the gap width 
may still appear surprising, but it can be under- 
stood with some physical intuition. In an inertial 
environment, it is pressure -and not viscosity - 
which makes a gas fill the void space. In a ro- 
tating disk the situation is different, because a 
radial pressure gradient simply adds or subtracts 
a force to the gravitational force exerted by the 
central star. This changes the angular velocity 
of rotation of the gas, without causing any radial 
transport of matter. Thus if the edges of the gap 
were circular, the pressure could not play any role 
in limiting the gap opening. However, as the gap 
edges are not circular, as shown in Fig. El the 
pressure gradient is not entirely in the radial di- 
rection, and thus it exerts a force with a non-null 
azimuthal component. This gives a net torque, 
and tends to fill the gap. 

4 Gap profiles 

In the last section, the pressure torque has been 
numerically computed in different cases. It has 
been shown that, when the disk is in equilib- 
rium, the pressure, gravity and viscous torques 
cancel out. This suggests that it should be pos- 
sible to compute a priori the shape of the gap by 
imposing that this equilibrium (jl2j) is respected. 
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Indeed, the viscous and pressure torques depend 
on the relative radial gradient of the azimuth- 
averaged density, whereas the gravity torque has 
no direct dependence on it. Therefore, on a given 
trajectory, there must be a value of this gradi- 
ent that corresponds to the exact equilibrium be- 
tween these three torques. 

Our semi-numerical algorithm for the computa- 
tion of this equilibrium value is described in ap- 
pendix B. The results are shown in Fig. [3 (crosses) 
and satisfactorily agree with the real values mea- 
sured in the corresponding numerical simulation 
(solid curve), i.e. the simulation from which the 
streamlines used by the algorithm have been ob- 
tained. 

The knowledge of the relative radial gradient 
of the azimuth-averaged density as a function of 
the radial distance enables us to construct a gap 
profile by simple numerical step by step integra- 
tion, starting from a boundary condition. In the 
secondary panel of Fig. [7| this integrated profile 
(dashed curve) is plotted against the real one from 
the considered simulation. The match between 
the two profiles is almost perfect, which again 
proves that the gap profile is set by the balance 
between the three torques due to gravity, viscos- 
ity, and pressure. 

4.1 An explicit equation for the 
gap profile 

We now wish to go beyond the semi-numerical 
algorithm of Appendix B and obtain an approxi- 
mate analytic expression for the pressure torque, 
to be used in an explicit differential equation for 
the gap profile. 

As we have seen above, for a given streamline, 
the absolute value of the pressure torque is an 
increasing function of the relative radial gradient 
of the azimuthally averaged surface density. Fur- 
thermore, in a disk with no density gradient, the 
pressure torque must be zero. Thus, we approx- 
imate the dependence of the pressure torque on 
the relative radial density gradient with a linear 



45 




r 



Figure 7: The crosses show the relative radial gra- 
dient of the surface density as computed by the al- 
gorithm described in Appendix B. The solid curve 
shows the same quantity, measured from the numer- 
ical simulation from which the streamlines used in 
the algorithm have been taken (aspect ratio = 0.05, 
viscosity = 10~^'^, planet mass = lO^"'). From the 
algorithm, the gap profile is computed (dashed curve 
in the little box), and compared to that obtained in 
the numerical simulation (solid curve). 

function : 

Before looking for a numerical approximation 
of the function a(r), we make two considerations 
on its functional dependence on the scale height 
of the disk and on the mass of the planet. 

First, because of Eq. (0), a(r) is necessarily pro- 
portional to Cg"^. As Cs is proportionnal to the 
scale height H, we can write a = {H/rya'{r). 

Second, in the limit of negligible viscosity, scal- 
ing the aspect ratio H/r proportionally to Rn/rp, 
and adopting Rh as basic unit of length, the 
equation of motion becomes independant of the 
planet mass (Korycansky and Papaloizou, 1996). 
Thus, if the disk aspect ratio scales with the 
planet Hill radius, the resulting surface density 
S at equilibrium is a function of A/Rh only. 
Consequently dS/(Sdr) is a function of A/Rh, 
divided by Rh- As the gravity torque tg is 
proportional to Rn'^r^A/ Rh)~'^ (see (II 1|) ), the 
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equilibrium tg = tp can hold if and only if 
a'(r) = rRn a"{A/RH) Vp^p^, where a"{A/RH) 
is a dimensionless function and the constant fac- 
tor Tpflp'^ stands for homogeneity reasons. 

To evaluate the function a", we use numerical 
simulations from which we measure the pressure 
torque and the relative radial gradient of the sur- 
face density. In practice, we consider two sim- 
ulations : (i) the reference one, with a Jupiter 
mass planet in a disk with aspect ratio = 0.05 
and viscosity = 10~^'^, which gives information 
for A/Rh in the range 2-7, and (ii) a similar 
simulation but with viscosity u = 10"^'^ which, 
because of its wider gap, allows us to better es- 
timate the asymptotic behavior of a" at large A. 
We find that a"{A/ Rh) can be approximately fit- 
ted by the function 
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Equation (fTT^ has been determined for the ex- 
ternal part of the disk (A > 0), outside of the 
horseshoe region. However, assuming that the 
streamlines are symmetric relative to the position 
of the planet, the same expression can be applied 
in the inner part of the disk, which justifies the 
absolute value of A. In fact, to represent the inner 
edge of the gap, just rotate Fig. Elby 180 degrees, 
and it becomes evident that a negative density 
gradient leads to a positive torque. 

Note that Eq. (fT!^ has been determined with 
reference to the streamlines corresponding to the 
case with q = 10"^, u = 10"^-^, and H/r = 5%. 
However, the exact shape of the streamlines de- 
pends on q, u, and H/r, even in rescaled coordi- 
nates. We neglect this dependence at this stage. 

Thus, we assume that (|T!?jl is valid for any value 
of u and H/r and a" depends on q only via Rh- 
This approximation has the advantage of provid- 
ing us an analytic expression for the computation 
of the gap profiles. In fact, the disk equilibrium 
equation (fT^ becomes : 



Rh dS 
S dr 
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Figure 8: Comparison of numerical results (plain 

lines) and analytic profiles given by Eq. 1)14(1 (dashed 

lines). Case 1 : reference case : q = 10~^, H/r = 5%, 

V = 10~^'^. Case 2: more viscous case: q = 10~^, 

H/r = 5%, u = 10"''. Case 3: less viscous case: 

q = 10-3, H/r = 5%, v = IQ'^-^ . Case 4: scaled 

case: q = SAO'^, H/r = 7.2%, u = W'^-^. Case 5: 

thicker disk case : q = 10"^, H/r = 10%, v = lO""^-^. 

Case 6: more massive case: q = 3.10"''^, H/r = 5%, 
= 10-5-5. 

with a" and tg given in Eq. (fT^ and (fTT|) respec- 
tively. 

The right hand side of Eq. ()14|) is independent 
of S and is an explicit function of r. This differen- 
tial equation can be integrated, once a boundary 
condition S(ro) is given. Unfortunately the inte- 
gral has no analytic solution, so that it has to be 
computed numerically. 

Figure |H1 shows comparisons of the gap profiles 
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obtained in numerical simulations with those ob- 
tained with the integration of Eq. (fT^ . for six 
different cases with different viscosities or aspect 
ratios and planetary mass (see figure caption for 
a list of parameters) . The comparisons are done 
only for the outer part of the disk, because in 
the inner part, the effect of the boundary condi- 
tion, not considered in our model, is too promi- 
nent in the numerical results. In the integration 
of Eq. (HD), S has been set equal to the value 
found in the numerical simulations at a point on 
the brink of the gap, so as to allow an easier com- 
parison between the numerical and semi-analytic 
density gradients at the edge of the gap. We re- 
mark that in case 1, the semi-analytic gap pro- 
file matches almost perfectly the numerical pro- 
file. This is not surprising, because this is the ref- 
erence case for which the streamlines have been 
computed, so that our expression (fT^ is virtually 
exact. 

In cases 2 and 3, we change the viscosity and 
keep the same planet mass and aspect ratio as in 
case 1. Now, the agreement between the numeri- 
cal and semi-analytic profiles is less good. In par- 
ticular, in the high- viscosity case, the real density 
gradient is shallower than the one we compute, 
while in the low-viscosity case it is steeper. This 
is because the real streamlines are not identical 
to those for which Eq. (fT3|l has been computed. 
In the more viscous case, the equilibrium in the 
disk is achieved with a weaker pressure torque. 
The distortion of the streamlines at the wake is 
dictated by the difference between the local grav- 
ity and pressure torques. Thus, a weaker pressure 
torque gives streamlines that are more distorted 
at the wake than in our reference case. But, as 
sketched in Fig. the more a streamline is dis- 
torted, the more efficient it is in producing a pres- 
sure torque from a radial surface density gradi- 
ent. Consequently the pressure torque required 
to set the equilibrium in the disk is achieved with 
a smaller density gradient than the one needed if 
the streamlines were as in the reference case. The 
opposite holds in the less viscous case. 

In case 4, we increase the mass of the planet and 



the disk aspect ratio, in a way such that H/Rh 
is the same as in case 1. The viscosity is also the 
same as in case 1, and the agreement between the 
model and the simulation is equally good. 

Finally, in case 5 and 6, we change H/Rh- In 
case 5, we keep the planet mass and viscosity 
of case 1 but increase the aspect ratio to 10% ; 
the model gap is quite narrower and shallower 
than the numerical one. In case 6, we use the 
same planet as in case 4 but with the aspect ra- 
tio and viscosity of case 1, which gives as smaller 
H/ Rh ; the gap that our model predicts is now 
slightly wider than the one obtained in the nu- 
merical calculation. The interpretation for the 
disagreements observed in cases 5 and 6 is the 
same as that offered for cases 2 and 3. 

4.2 Note on disk evolution during 
gap opening 

Figure ISl shows significant differences in the value 
of S between the numerical simulation and the 
analytic expression. However, we stress that a 
large difference in S can correspond to almost no 
difference in the relative slope (^^)- For in- 
stance, in the case with viscosity equal to 10~^ 
(case 2), the surface density profiles for r > 1.6 
seem quite different, but in fact, they have the 
same relative slope. 

As we have shown above, it is the relative slope 
(^^) that sets the equilibrium. This equilib- 
rium must be reached quickly, on a time scale 
independent of viscosity. In fact, in absence of 
equilibrium, the fluid elements are displaced ra- 
dially over a synodic period and the trajectories 
are not periodic. This corresponds to the open- 
ing of the gap. Then, once the equilibrium is 
almost set, the value of S can still significantly 
evolve on a long (viscous) time scale, but keeping 
(^^) essentially unchanged. This fact is illus- 
trated on Fig. El which compares the evolution of 
(^^) (top panel) with the evolution of S (bot- 
tom panel) for a weakly viscous case [v = 10~^'^). 
This behavior explains why, when simulating the 
gap opening in low viscosity disks, the surface 
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Figure 9: Evolution with time of the profile of the 
external edge of the gap opened by a Jupiter mass 
planet in a 5% aspect ratio disk, with a 10~^'^ vis- 
cosity. Top panel shows the relative radial gradient 
of the density ; the fact that the curves overlap ar- 
gues that the equilibrium function has been rapidly 
reached. Bottom panel shows the evolution of the 
surface density at the same times ; a 'bump' appears, 
which is a consequence of the matter removed from 
the gap, and it is eroded on a viscous time scale. This 
happens without modifying substantially the relative 
slope. 



density profile seems to have attained a station- 
ary solution within a limited number of plane- 
tary orbits, despite that the presence a prominent 
'bump' at the outer edge of the gap indicates that 
there is still room for evolution (see Fig. Q). 



In the inner disk, once the gap profile is set in 
terms of relative slope, we expect that the density 
E decreases on a viscous timescale, because of the 
accretion on the central star. In the approxima- 
tion of a fixed planet, this viscous evolution would 
lead to the formation of an inner hole in the disk, 
extended up to the planet position. 



5 Dependence of gap pro- 
files on viscosity and as- 
pect ratio 

In the previous section, we have presented a semi- 
analytic method to compute gap profiles. It gave 
overall satisfactory results, as shown in Fig. |H| In 
this section, we use our method to explore the 
dependence of the gap profile on the two key pa- 
rameters of the disk : viscosity v and aspect ratio 
Hjr. In particular, we wish to revisit, with a 
unitary approach, the gap opening criteria men- 
tioned in section |21: 

(i) the viscosity needs to be smaller than a thresh- 
old value. According to Eq. in our case of a 
Jupiter mass planet this value is z/crit ~ 10~^. 

(ii) The disk height at the location of the planet 
needs to be smaller than ~ Rn. For a Jupiter 
mass planet it corresponds to an aspect ratio ~ 
7%. 

In the computation of the gap profiles by in- 
tegration of Eq. (|T^ . two problems are encoun- 
tered. 

First, a boundary condition is needed. This 
choice is arbitrary, but in principle it should be 
consistent with the steady state of the disk. How- 
ever, the steady state is an academic concept 
which exists only if the density is kept fixed some- 
where in the disk, otherwise the disk spreads to 
infinity following Lynden-Bell & Pringle (1974) 
equation. In our numerical code, the surface den- 
sity is kept equal to the unperturbed value at the 
outer boundary of the grid (r = 3). Thus, for the 
solutions of Eq. (fT^ presented in Figs. [TUl and 
ITU we impose S(r = 3) = l/-\/3- In this way, 
our solution should correspond to the steady state 
solution that the code would converge to. More- 
over, this choice allows a direct comparison of our 
profiles with those obtained with Varniere et al. 
model, illustrated in the top panel of Fig. [T]using 
the same boundary condition. 

The steady state solution provided by the nu- 
merical simulation does not depend on the size of 
the grid, provided that the boundaries are suffi- 
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ciently far from the planet (negligible differential 
planetary and pressure torques, or equivalently, 
negligible wave carried angular momentum flux, 
as in our nominal case - see Appendix C ). This 
required size increases with decreasing viscosity 
because the radial range over which the wave is 
damped increases. 

Thus, in a very low viscosity case, the steady 
state solution obtained by the numerical simula- 
tion over an extended disk would be different from 
the model profiles given on Fig. ^| However, our 
model profiles would bound the gap observed in 
the numerical simulation as long as the normal- 
ized surface density at r = 3 does not decrease 
below 1/a/3. In such low viscosity cases, this 
happens after an exceedingly long time. Thus, 
we claim that our model profiles are significant 
for the description of gaps in realistic disks. 

The second problem concerns the treatment of 
the horseshoe region. The gravity and pressure 
torques, tg and tp, are considered null in the 
horseshoe region |A| < 2Rh- The depth of the 
gap is thus set by the value of the density at 
rp + 2Rh. At the edges of the gap, the slope 
is very steep, so that a little change in the as- 
sumed width of the horseshoe region leads to a 
major change in the gap depth. This is a limi- 
tation of our results from a quantitative point of 
view. Though, it does not change the qualitative 
evolution of the gap profiles with respect to the 
disk parameters. 

This sensitivity to the width of the horseshoe 
region is also a problem for the construction of the 
surface density profile in the inner disk. The in- 
tegration for the inner disk starts from rp — 2 Rfj 
down to r = 0, with the density at the bottom 
of the gap acting as the boundary condition. In 
principle, if the gap profile is symmetric, the er- 
rors at the right hand side and left hand side bor- 
ders of the gap compensate each other : the value 
of the surface density at the bottom of the gap is 
not quantitatively correct, but the density profile 
in the inner disk is realistic. 

Figure ITUl shows the results of our semi-analytic 
calculation for a fixed value of the aspect ratio 





Aspect Ratio = 0.05 




log(nu) = -7.5 




log(nu) = -6.5 - 




log(nu) = -6.0 




log(nu) = -5.5 




log(nu) = -5.0 - 




log(nu) = -4.5 




log(nu) — -4.0 




\ log(nu) = -3.14 . 




C\ sqrt ( 1 / X ) 







0.5 1 1.5 2 2.5 3 

r 

Figure 10: Analytical gap profiles given by Eq. ([TUl 
for different viscosities. The gap deepens as viscosity 
decreases, but its width remains bounded, even for 
u = 0. 
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Figure 11: Analytical gap profiles given by Eq. (fTH) 
for different aspect ratios. The gap deepens as the 
aspect ratio decreases. 

(5%) and different values of the viscosity u (from 
to 10~^ in normalized units). Figure fTTl keeps 
the viscosity u = 10~^'^ and explores the depen- 
dence of the gap profile on the disk aspect ratio 
(from to 30%). The plotted curves naturally 
order themselves from bottom to top, from the 
less viscous case (respectively the smallest aspect 
ratio) to the most viscous case (respectively the 
biggest aspect ratio). Notice that this progres- 
sion does not represent an evolution with time, 
but different steady state gap profiles, for differ- 
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ent parameters. 

dependence on viscosity : First of all, we 
remark on Fig. El that the different shapes of 
the gaps qualitatively agree with those computed 
with numerical simulations, shown in the bot- 
tom panel of Fig. [T] Indeed, not only do we get 
deeper and wider gaps as viscosity decreases, but 
we also correctly reproduce the limited gap width 
achieved in the small viscosity cases. This means 
we have solved the problem that initially moti- 
vated our investigation. 

As viscosity increases, the gap is filled with gas, 
and the profile tends to the unperturbed profile 
set by the sole viscous effects: S oc I/a/t (see 
section 121). Nevertheless, it is hard to determine 
a precise threshold value for gap opening, for at 
least two reasons. The first one is that the gap 
profiles have a smooth dependence on the viscos- 
ity. The concept of threshold viscosity for gap 
opening does not hold. The gap gradually in- 
creases in depth over a range of viscosities of 
about one order of magnitude. The second rea- 
son is that the depth of our gaps is very sensitive 
to the assumed width of the horseshoe region, as 
discussed above. Thus, there is some uncertainty 
on the value of the viscosity that makes the gap 
become only a dip. Assuming that a gap is 
opened if the surface density falls below 10% of 
the unperturbed density, we find that z/crit ~ 10"^ 
We remind that the 'classical' threshold for gap 
opening is Ucrit ~ 10~^. However, the numerical 
experiments in Fig. [T] (bottom panel) also suggest 
that S ~ 0.1 at the bottom of the gap is achieved 
for u ~ 10~^. 

dependence on aspect ratio : Consider now 
Fig. ^2 We see a smooth evolution from deep 
gaps to shallow or inexistent gaps with increasing 
aspect ratios. This is easy to understand, be- 
cause [H/rY is a multiplicative coefficient in the 
expression of the pressure torque (see section H}. 
Therefore, the larger (H/r), the shallower needs 
to be the relative slope at the edge of the gap to 
achieve the equilibrium ()12|1 . As in the previous 



case, it is not possible to determine a threshold 
value of {H/r) for gap opening, but we find that 
the 'classical' value {H/r)crit ~ 0.07 corresponds 
to about 90% depletion in the gap. 

More generally, with our approach we find that 
the viscosity required to fill the gap is a decreasing 
function of the aspect ratio. If the aspect ratio 
is too large, the gap cannot be opened whatever 
the viscosity. To our knowledge, this is the first 
time that an analytic approach gives the correct 
description of the evolution of the gap profile with 
respect to both disk viscosity and aspect ratio. 

6 A new generalized crite- 
rion for gap opening 

To go beyond the qualitative considerations of 
the previous section, we try to generalize the gap 
opening criterion with an expression that involves 
simultaneously the three main parameters of the 
problem: mass of the planet, scale height of the 
disk and viscosity. 

We start with a few considerations on two lim- 
iting cases. In the zero viscosity limit, as we have 
seen in section IH changing the scale height of 
the disk in proportion to the Hill radius of the 
planet preserves the gap profile in scaled units 
A/Rh. This means that, whatever depth thresh- 
old is adopted for the definition of 'gap', the 
threshold value of H for gap opening in the zero 
viscosity limit, Hq, is proportional to Rh '■ 

1 

Hq oc Rh oc gs . 

In the infinitely thin disk limit {H/r — > 0), the 
disk equilibrium is set by the equation tg = ty. At 
the border of the gap where the slope of the sur- 
face density is relevant, ty is proportional to z^^^^. 
By changing the mass of the planet, the gravity 
torque changes proportionally to RH^r{A./ Rh)~^ ■ 
If the viscosity v is changed proportionally to 
Rh^ oc g, then the surface density profile S re- 
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mains an invariant function of A/Rh- Thus, in- 
dependently of the adopted definition of 'gap' as 
in the previous case, the threshold viscosity for 
gap opening in the infinitely thin disk limit, z/q, 
scales proportionally to q. This is consistent with 
the gap opening criterion given in Bryden et al. 
(1999). 

I/O oc g 

We now come to the general case where neither 
H nor V are null. From the considerations above 
and Eq. (fT^ it is evident that a change in the 
planet mass q can give an invariant surface den- 
sity profile in scaled units A/Rh provided that H 
is changed proportionally to and v is changed 
proportionally to q. 

The most complicated case that remains to be 
analyzed is that where q is constant, but H and v 
are changed. It is evident from Eq. (fT^ that it is 
not possible to have an invariant surface density 
profile by decreasing H and increasing v or vice- 
versa. The question is therefore how to keep the 
central gap depth invariant, despite changes in 
the gap profile. We answer this question using our 
semi-analytic calculation of gap profiles, based on 
the integration of Eq. ()14|) . For this, we define - 
arbitrarily - that the minimal depth that defines a 
gap is 1/10 of the unperturbed disk density at r = 
rp. Figure IT^ shows as bold lines, for six different 
values of the planet mass, the relationships H vs. 
V that preserve such central gap depth. 

As one can see, these relationships are almost 
linear. 

We can fit each one with a relation of type 
H/Hq + zz/z/q = 1, where Hq and vq have been 
defined above. As we have z/q oc g and Hq oc g^/'^, 
we can derive a general relation involving H, v 
and g that approximately describes all the curves 
plotted on Fig. ^[ and thus a general criterion 
for gap opening. Denoting by Tl the Reynolds 
number rpVLp/p, we find that a gap is opened if 
g, H and IZ satisfy the following inequality : 
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Figure 12: The bold curves represent the values of 
(H/r) that make the gap depth be 10% of the unper- 
turbed density for given values of u. They have been 
computed from the solutions of Eq. (|14|) . Each curve 
corresponds to a planet mass. The thin lines repre- 
sent our linear approximations given by Eq. (|15jl for 
the corresponding planet mass. 
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7 Conclusion 

In this paper, we have analyzed in detail the pro- 
cess of gap opening in proto-planetary disks. In 
this respect, a key problem is to calculate which 
fraction of the torque exerted by the planet is 
locally deposited in the disk and which fraction 
is transported away by pressure supported waves. 
We have shown that the angular momentum evac- 
uated by the waves can be computed as a pres- 
sure torque. We found that the steady state of 
the disk is set by the equilibrium among the total 
gravity torque, the viscous torque, and the pres- 
sure torque. From this consideration, we have 
built a semi- analytical algorithm that, given vis- 
cosity and aspect ratio, provides the equilibrium 
profile of the surface density of the disk, enabling 
us to explore the gap shape for a large range of 
parameters. 

This work has two types of application. It can 
be used to achieve a first realistic estimate of the 
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width and depth of gaps in various situations, in 
view of the future high resolution observations of 
proto-planetary disks (with ALMA or the SKA 
projects). It can also give equilibrium gap pro- 
files to be used starting condition in numer- 
ical simulations if one wants to avoid the inter- 
mediate, cpu-consuming phase which leads to the 
steady state. 

Our work is not fully analytic. Indeed our final 
equation (fT^ involves a function a" which we ap- 
proximated by the ansatz function (fT^ . with co- 
efficients determined with respect to a reference 
numerical simulation. Also the gravity torque 
flll|) has been refined using fits to the reference 
numerical results. As a consequence, if our model 
matches the results of the reference numerical 
simulation, it still is in satisfactory agreement 
with the results of other numerical simulations. 

Moreover, the equilibrium profile that we ob- 
tain corresponds to the equilibrium configuration 
of the disk at infinite time in presence of a non- 
migrating planet, which is evidently an ideal case. 
Our model is two-dimensional, intended to ap- 
proximate the behavior of a vertically isothermal 
3D disk; in a 3D, thermally stratified disk, the 
density waves would not propagate exactly the 
same way (Bate et ai, 2003) and consequently the 
pressure torque is expected to be somewhat dif- 
ferent. Finally, we have assumed a constant kine- 
matic viscosity ; in reality, in the regions where 
the perturbations are nonlinear, the effective vis- 
cosity depends on the local planet's gravitational 
torque (Goodman and Rafikov, 2001), although 
this dependance may be weak (Papaloizou et ai, 
2004). 

In spite of these limitations, our work clearly 
demonstrates the fundamental role of the pressure 
in setting the equilibrium of the disk. Moreover, 
it gives a correct, nearly quantitative, description 
of the evolution of the gap profile with respect to 
the key parameters of the problem : planet mass, 
viscosity and aspect ratio. From this we derive a 
new general criterion for gap opening, involving 
simultaneously these three parameters. 

Our work shows why the width of the gap is 



bounded even in the case with very small viscos- 
ity, which was the open problem that originally 
motivated our work. It provides a conceptual uni- 
fication of the two classically, but independently 
derived, criteria for gap opening, based on thresh- 
old viscosity and aspect ratio. 
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8 Appendix 

A: trajectories and streamlines 
computation. 

In a steady state, trajectories and streamlines co- 
incide. Computing the streamlines is then equiv- 
alent to computing the trajectories. In principle, 
to compute a trajectory it is enough to integrate 
the velocity field. The latter is defined on the 
grid and output by the code, from which the ve- 
locity at any point of the disk can be computed 
by interpolation. However, using this procedure, 
the resulting trajectories would in general not be 
periodic, as a consequence of the accumulation of 
the integration and interpolation errors. This is 
a serious problem, because the loss of periodic- 
ity would introduce a spurious change of angular 
momentum, namely a spurious torque. 

To obtain perfectly periodic streamlines, we 
used the following algorithm, that for simplicity 
we detail for the outer part of the disk (r > rp). 
We first compute a trajectory from (r = Tq, 6' = 
vr) to ^ = —7T by simple numerical integration of 
the velocity field (the integration runs from vr to 
—71 because tq > rp, so that the fiuid element 
rotates clockwise in the corotating frame). This 
gives a first curve r^^\9), defined on the interval 
[— 7r,7r]. By definition A^^tc) = tq, but r^^^{—7i) 
is in general close but not equal to tq, because 
of numerical errors, as said above. On this tra- 
jectory, we calculate ^{r^^\e),e) = f'^^\e). This 

is a pseudo-derivative of r«, i. e. the slope of 
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the tangent to the curve according to the velocity 
field. It should be equal to dA^^dO, but is not 
exactly equal because of the numerical errors in 
the computation of r^^\ Then, we compute the 
Fourier coefficients fn^ of f''^\6). The first one 
/o^'* is real, and corresponds to the mean of f^^\ 
namely to a radial drift. It is not zero as r^^^ is 
not exactly periodic, and therefore we set it to 
zero. The pseudo-derivative of r^^^ with respect 
to 9 is thus modified. To get back to a trajec- 
tory, we integrate this modified pseudo-derivative. 
We denote the new trajectory by r^'^\6). This 
trajectory is periodic by construction as its ze- 
roth order Fourier coefficient is null. From r^'^\ 
we repeat the algorithm to find r*-^\ and so on, 
until the algorithm converges to a fixed point. 
This fixed point r{6) is a periodic trajectory by 
construction. It fits the velocity field, as it veri- 
fies = ^{r{6),6), provided that the zeroth 
order Fourier coefficient of its pseudo-derivative 
is negligible. If it isn't, it means that the algo- 
rithm failed. This happens in particular if the 
real streamlines are not periodic because a steady 
state has not been reached yet. 

In practice, for the implementation of this algo- 
rithm, we used simulations computed over a grid 
with a larger resolution than that used in sec- 
tion El More precisely, we have used 512 cells 
in radius and 1024 in azimuth. The number of 
points used to compute the Fourier coefficients 
of the pseudo-derivative was 1024 too. In all 
cases, the algorithm explained above converged, 
and the zeroth order Fourier coefficient of the final 
pseudo-derivative was negligible (less than 10~^ in 
our normalized units, even lO"'^ for all trajecto- 
ries with r(7r) > 1.2). 

B: a semi- numerical algorithm for 
the calculation of the equilibrium 
surface density slope. 

We present an algorithm that, given the shape 
of the streamlines, computes the relative surface 
density radial gradient that ensures the equilib- 



rium condition (fT^ . This is done in two steps. 
First, we design a procedure that evaluates tp on 
each streamline, for any given value of Sec- 
ond, we solve numerically the implicit equation 
for given by Eq. (IT^ . 

First step : computation of the torques 

The streamlines are ordered with respect to in- 
creasing distance to the central star, so that 
ri{9) < rj+i(^) for every i, 9. We call the ith 
streamtube the zone around the zth streamline : 
{ (rj_i +ri)/2 <r < {ri + ri+i)/2 }. A total mass 
rrii or mean density Sj can be imposed to be car- 
ried by a given streamtube i. Because the steady 
state is reached, the fiux of matter in streamtube 
i is constant with respect to time and azimuth, 
and is Fj = rrii/Ti, where Tj is the synodic period 
along the streamline. Thus, the mass has to be 
distributed in the streamtube in such a way that 
the fiux 

F{9) = J:ir,{9),9)ve{r0),9)[r,+,i9) - r,_,i9)]/2 

(16) 

is equal to Fi for all 9. The azimuthal speed 
Ve{fi{9)i 9) can be obtained by interpolation from 
the output of the numerical code ; the local den- 
sity Tj{ri{9),9) is therefore the only unknown in 
Eq. (fTB|) , so that one has : 

nr^{^). 9) = 2Fi/v0{r,{9), 9) [r,+^{9) - r,_^{9)] . 

(17) 

Any relative radial density gradient 
(l/S)(dS/dr) around the ith streamline can be 
created by imposing appropriate values for Sj+i 
and The masses rrii+i and carried by 

the streamlines are obtained by multiplying the 
mean surface densities by the areas of the stream 
tubes. Then, the local densities are computed 
using Eq. (fT7|). 

Once the streamlines and the local densities are 
known, the numerical computation of the pres- 
sure torque can be done using Eq. (fTUj) . The 
partial derivative of the density with respect to 
the azimuth is delicate to compute. Indeed, from 
Eq. (fTTjl we know S(rj(6'), 9) only on a discrete set 
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of values ri{9). To compute {dT,/d9) at the loca- 
tion {ri{6),6) we need to know J]{ri{6),6 ± 66), 
for some small 66. This is computed by in- 
terpolation between Tj{rj{6 ± 66), 6 ± 66) and 
I^{rj+i{6±66), 6 ±66), where the jth streamline is 
chosen such that rj{6±66) < ri{6) < rj+i{6±66). 

The viscous and gravity torques are given by 
expressions (j3)) and (lllj) respectively, with r = 
r(7r). Expression Q is preferred to expression 
Q, despite of its limitations in the very vicinity 
of the planet (see Fig. ^ because it is simple and 
explicit. 



Second step : computation of the gap pro- 
file To obtain the density profile we impose that 
the sum of the three torques vanishes on every 
streamline. Thus, for each streamline, the goal 
is to find the value of the relative surface den- 
sity slope s = (^^) that makes the total torque 
^totai(s) = (tp + tu + tg) equal to zero. As the 
pressure torque is numerically computed, the so- 
lution can be found only numerically. We use a 
secant method algorithm, described next. 

A first value Sq of s is arbitrarily chosen (typ- 
ically or the solution found on a neighboring 
streamline). Then, another value si is taken (for 
instance So-t-lOOttotaK-So))- The secant method al- 
gorithm is then used. The value chosen for S2 is : 

Sl-itotal(Sl)-(Sl-So)/ [itotal(Sl) - ttotal(so)]- A Se- 

quence (s„)„=o,i,... is build this way. At each step, 
the tested value is : Sn = s„_i-ttotai(s„-i).(s„__i- 
Sn-2)/ [^totai(sn-i) -^totai(s„--2)]- The scqucuce 

converges to SequH, such that ttotal(Sequil) = 0. 

We stop when we reach a value for s that makes 
|^totai(s)| smaller than 10~Hg, and we take that as 
our solution for (^^)equii- 

With this procedure, we get (^^)equii for each 
streamline or, equivalently, each rj(7r). It repre- 
sents a data point for the relative radial derivative 
of the density, shown as a cross on Fig. [7| 



C: Flux of angular momen- 
tum. 

The flux of angular momentum has to be evalu- 
ated in a frame in which angular momentum is 
conserved. This is not the case for the frame 
centered on the primary (which is accelerated), 
whereas it is the case in the non-rotating frame 
centered on the barycenter G of the system (star 
plus planet plus disk), which is inertial. One 
therefore needs to evaluate the following expres- 
sion : 

Fh= / {T.rv'e)v'^rd6, (18) 
Jo 

where v'q and v'^ are the perturbed azimuthal 
and radial velocities in the G centered frame : 
v'q = V0—V0 and v'^, = v^—Vr, the barred quantities 
denoting the averages over the circle of integra- 
tion. 

We assume that g <^ 1. We remark that the 
perturbed quantities are proportionnal to q, and 
Fh to q^. Then, to compute (fTHj) from the veloc- 
ities output by the code, we need a sequence of 
transformations. Neglecting terms that will give 
corrections of order q^ in F^, this reduces to two 
transformations on Vr and Ve : 

(i) The velocity of G in the heliocentric frame 
has to be substracted. In polar coordinates 
centred on the star, it is to first order in q : 
v{G) = qrpilp[sm{6 — 6p),cos{6 — 6p)), where 
the subscript p refers to the planet. 

(ii) The radial and azimuthal components of a 
fiuid element velocity are different in the he- 
liocentric and barycentric frames. For any 
vector X = (Xr,Xe) in the heliocentric 
frame, the radial component in the barycen- 
tric frame is written, to first order in q, as : 
Xr — Xgq{rp/r) sm{6 — 6p). Similarly, the az- 
imuthal component of X in the barycentric 
frame is: X0 + Xrq{rp/r) sm{6 — 6p). We 
stress that the radial component of the ve- 
locity of a fiuid element is proportional to q, 
so that the above correction on the azimuthal 
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component is second order in q and will be 
neglected. 

The application of (i) and (ii) give the following 
expression for Fh '■ 

Fh = [r {v'g - qrpVLp cos{e - 9p)) x 

Jo 

- g^(rfip + Ve) sm{e - 0^)) ] rdO , (19) 

where S is the mean density on the circle and all 
the quantities are the ones output by the code in 
the heliocentric frame. This corresponds to the 
flux of angular momentum through the circle of 
radius r, due exclusively to the wave launched by 
the planet. 

The assumption g <^ 1 has allowed us to ne- 
glect the following effects : 

(a) The density S should be evaluated along the 
circle, but E = S + E', and S' oc gE -C S, so 
that S can be replaced by S in the integral. 

(b) The circle of radius r centered on G differs 
from the circle of radius r centered on the 
star. As the distance between the two cir- 
cles is proportional to q, this only introduces 
negligible modifications in the value of every 
quantity. 

(c) In the previous calculations, G corresponds 
to the barycenter of the star-planet system, 
and not of the whole system including the 
disk. The latter is initially axisymmetric, 
and the perturbations are proportional to q. 
As the mass of the disk is also of the order of 
the mass of the planet, the influence of the 
disk on the barycenter position is negligible. 

We computed the flux Ffj on our reference sim- 
ulation {q = 10-^ u = 10-^•^ H/r = 0.05) us- 
ing Eq. p9|) . In Fig. 1131 Fh{t) is plotted as a 
bold plain line, whereas the total gravity torque 
Tg computed on the annulus between the planet 
orbit and the circle of radius r is shown clS Qb bold 
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Figure 13: Angular momentum flux carried by the 
wave launched by the planet (bold plain line, corre- 
sponding to Eq. (dJ), compared to the total gravity 
torque (bold dashed line) as functions of the distance 
to central star. The difference is plotted thin 
dot-dashed line. 



dashed line. The gravity torque is computed us- 
ing the direct terms due to the planet {GMp/cP, 
d being the distance between the planet and the 
considered point) and to the star (GM^/r^), as it 
is evaluated in an inertial frame. The difference 
between Fh and Tg is the thin dot-dashed hue ; it 
corresponds to the cumulative locally deposited 
gravity torque (i. e. the fraction of the gravity 
torque that is not evacuated by the pressure sup- 
ported wave). The wave carries an increasing flux 
near the planet (in the zone {1.15 ^ r < 1-5}), 
and takes away a large fraction of the gravity 
torque ; this corresponds to the raduis where the 
pressure torque tp appears to be of the same order 
as the gravity torque tg (see Fig. This angu- 
lar momentum is then deposited further from the 
planet, in particular in the {1.5 ^ r < 2} region, 
where Fh{t) sharply decreases. Beyond r ~ 2 the 
flux vanishes. At the outer boundary of our grid, 
the flux of angular momentum taken away by the 
wave is negligible with respect to the total gravity 
torque. 

This shows that the outer boundary of the grid 
is sufficiently far from the planet so that the angu- 
lar momentum transfer from the wake to the disk 
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is correctly described, while the angular momen- 
tum leakage at the outer boundary is negligible. 
Thus, we conclude that our simulations arc real- 
istic, and our gap profiles correspond to steady 
states in the non-migrating planet hypothesis. 
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